#####################################################
#####################################################
## Seminář o práci v R Studiu - Část čtvrtá
## 
## Regrese
#####################################################
#####################################################

# Program:
# 1. Import a úprava dat
# 2. Jednoduchá regresní analýza
# 3. Vícenásobná regresní analýza
# 4. Export výsledků regresní analýzy
# 5. Reakce na porušení předpokladů (transformace dat)





###################################################
## 1. Import a úprava dat
###################################################

# Analýza volební účasti ve volbách do Evropského parlamentu 2019:
eurovolby <- read.csv("data/eurovolby_2019.csv")
# zdroj: https://volby.cz/opendata/ep2019/ep2019_opendata.htm 
eurovolby <- eurovolby[eurovolby$TYP_OBCE != "Městská část",] # vyřadíme městské části

# Počet obyvatel v obcích k 1. lednu 2019:
library(readxl)
obyvatel <- read_excel("data/obyvatele_0119.xlsx", sheet = 1)
# zdroj: https://www.czso.cz/csu/czso/pocet-obyvatel-v-obcich-za0wri436p

# Nezaměstnanost v obcích v dubnu 2019:
nezamestnanost <- read_excel("data/nezamestnanost_0419.xlsx", sheet = 1)
# zdroj: MPSV


# Kombinace datových setů:
eurovolby <- merge(x = eurovolby, y = nezamestnanost[,c("KOD","NEZAMESTNANYCH")], by.x = "KODZASTUP", by.y = "KOD")
eurovolby <- merge(x = eurovolby, y = obyvatel[,c("KODZASTUP","OBYVATEL","ZENY","VEK")], by = "KODZASTUP")

# Úprava proměnné:
eurovolby$ZENY_PROC <- eurovolby$ZENY/eurovolby$OBYVATEL*100



###################################################
## 2. Jednoduchá regresní analýza
###################################################

# Naší alternativní hypotézou je, že míra nezaměstnanosti v obcích měla vliv na volební účast.
# Co je nulová hypotéza?

# Nejprve se podívejme na rozložení dat:
hist(eurovolby$NEZAMESTNANYCH, breaks = 50)
hist(eurovolby$UCAST, breaks = 50)

plot(eurovolby$NEZAMESTNANYCH, eurovolby$UCAST) # kontrola vztahu dat

# Regresní analýza:
reg_ucast_m1 <- lm(formula = UCAST ~ NEZAMESTNANYCH, data = eurovolby)
summary(reg_ucast_m1)

# Ve výsledku vidíme hodnotu průniku (32.60) a koeficientu nezaměstnanosti (-0.86).
# Obě hodnoty mají určitou chybu a t-hodnoty.
# Zásadní je poslední sloupec Pr(>|t|) - u koeficientu vidíme, že je p-hodnota v podstatě nulová.
# Můžeme tedy téměř se 100% pravděpodobností zamítnout nulovou hypotézu.
# Dále vidíme vysvětlovací schopnost modelu, která činí 3.2 % vsysvětleného rozptylu závisle proměnné.
# F-statistický test ukazuje na celkovou kvalitu modelu - čím je výsledek vyšší, tím lépe.

# Zakreslení výsledku:
plot(eurovolby$NEZAMESTNANYCH, eurovolby$UCAST,
     xlab = "Míra nezaměstnanosti v obci (%)", ylab = "Volební účast (%)")
abline(reg_ucast_m1, col = "red", lwd = 2)

# Intervaly spolehlivosti:
confint(reg_ucast_m1, level = 0.99)

# Zakreslení:
abline(confint(reg_ucast_m1, level = 0.99)[,1], lt = "dashed", col = "red")
abline(confint(reg_ucast_m1, level = 0.99)[,2], lt = "dashed", col = "red")

abline(confint(reg_ucast_m1, level = 0.95)[,1], lt = "dashed", col = "darkgreen") # menší interval
abline(confint(reg_ucast_m1, level = 0.95)[,2], lt = "dashed", col = "darkgreen")

legend("topright", legend=c("99% interval spolehlivosti", "95% interval spolehlivosti"), # legenda
       col=c("red", "darkgreen"), lty="dashed", cex=0.8)


# Všimněte si, že regresní přímka protíná křížení průměrů obou proměnných:
abline(v = mean(eurovolby$NEZAMESTNANYCH), col = "yellow", lwd = 2)
abline(h = mean(eurovolby$UCAST), col = "yellow", lwd = 2)

# Proč tomu tak je?



###################################################
##  Cvičení I: Přišli k volbám spíše starší voliči?
#   a) Pomocí regresní analýzy zjistěte vztah mezi věkem a volební účastí.
#   b) Je vztah statisticky významný? Pokud ano, popište ho...





reg_ucast_vek_m1 <- lm(formula = UCAST ~ VEK, data = eurovolby)
summary(reg_ucast_vek_m1)



###################################################
## 3. Vícenásobná regresní analýza
###################################################

# Naší alternativní hypotézou dále je, že nezaměstnanost není jediným faktorem, 
# který ovlivnil volební účast ve volbách do Evropského parlamentu.
# Dalšími faktory, u kterých předpokládáme vliv, jsou:
# velikost obce, průměrný věk v obci a procentuální zastoupení žen.

# Znovu se nejprve podívejme na rozložení dat:
hist(eurovolby$UCAST, breaks = 50)
hist(eurovolby$NEZAMESTNANYCH, breaks = 50)
hist(eurovolby$VEK, breaks = 50)
hist(eurovolby$ZENY_PROC, breaks = 50)
hist(eurovolby$OBYVATEL, breaks = 50)


# Regresní analýza:
reg_ucast_m2 <- lm(formula = UCAST ~ NEZAMESTNANYCH
                   + OBYVATEL + VEK + ZENY_PROC, data = eurovolby)
summary(reg_ucast_m2)

# Jaké koeficienty nezávisle proměnných jsou statisticky signifikantní na úrovni 99 %?

# Podíváme se na splnění předpokladů:
install.packages("car")
library(car)
plot(eurovolby$NEZAMESTNANYCH, eurovolby$UCAST) # lineární vztah proměnných
vif(reg_ucast_m2) # test multikolinearity
cooks.distance(reg_ucast_m2)[cooks.distance(reg_ucast_m2)>1] # test odlehlých hodnot
par(mfrow=c(2,2))
plot(reg_ucast_m2) # test homoskedasticity a normální distribuce reziduí
dev.off()
durbinWatsonTest(reg_ucast_m2) # test nezávislosti reziduí


# Simulace intervalů spolehlivosti:
range(eurovolby$NEZAMESTNANYCH) # nezaměstnanost se pohybuje cca mezi 0 a 20 %
# Nejprve vytvoříme simulační matici nezávisle a kontrolních proměnných.
# Kontrolní proměnné nastavujeme do pozice středních či co nejvíce frekventovaných hodnot.
# Jde tedy o simulaci co nejběžnějšího a tudíž nejpravděpodobnějšího scénáře:
ci <- data.frame(NEZAMESTNANYCH = seq(0,20,1), OBYVATEL = mean(eurovolby$OBYVATEL),
                 VEK = mean(eurovolby$VEK), ZENY_PROC = mean(eurovolby$ZENY_PROC))
# Následně odhalíme hodnoty závisle proměnné
sim <- predict(reg_ucast_m2, newdata = ci, interval = "confidence", level = 0.99)

# Výsledek zakreslíme:
plot(eurovolby$NEZAMESTNANYCH, eurovolby$UCAST,
     ylab = "Volební účast (%)", xlab = "Míra nezaměstnanosti (%)",
     main = "Vztah míry nezaměstnanosti a volební účasti (komplexní model)")
lines(c(0:20),sim[,1], lt = "solid", lwd = 2, col = "blue")
polygon(c(c(0:20),rev(c(0:20))), c(sim[,2],rev(sim[,3])), 
        col = adjustcolor("blue",alpha = 0.2),  border = FALSE) # zakreslení polygonu
lines(c(0:20),sim[,2], lt = "dashed", lwd = 1, col = "blue")
lines(c(0:20),sim[,3], lt = "dashed", lwd = 1, col = "blue")
abline(reg_ucast_m1, col = "red", lwd = 2) # srovnání s jednoduchou regresí



###################################################
##  Cvičení II: A co moje obec?
#   a) Jak se vyvine volební účast v mojí obci, když stoupne z 5 na 9 %?
#      U nás žije 673 obyvatel, kterým je průměrně 37 a z nich je 51 % žen.





ci <- data.frame(NEZAMESTNANYCH = c(5,9), OBYVATEL = 673,
                 VEK = 37, ZENY_PROC = 51)
predict(reg_ucast_m2, newdata = ci, interval = "confidence", level = 0.99)



###################################################
## 4. Export výsledků regresní analýzy
###################################################

# Existuje několik možností - jednou z nejlepších je funkce stargazer.
install.packages("stargazer")
library(stargazer)
help(stargazer)

# Regresní tabulku zapisuje například jako formát html.
# Následně stačí stejný soubor otevřít ve wordu a využít příslušnou tabulku.
stargazer(reg_ucast_m1, type = "html", out = "export/reg_ucast.html")

# Exportovat je možné i několik modelů současně.
stargazer(list(reg_ucast_m1, reg_ucast_m2), type = "html", out = "export/reg_ucast.html")



###################################################
## 5. Reakce na porušení předpokladů (transformace dat)
###################################################

# Podívejme se znovu na distribuci dat u počtu obyvatel v obcích.
hist(eurovolby$OBYVATEL[eurovolby$OBYVATEL < 10000], breaks = 50) # data jsou pozitivně zešikmená 

# Předpokladem lineární regresní analýzy ale je, že distribuce dat se co nejvíce blíží normálnímu rozložení.
# Je proto třeba přistoupit k transformaci dat -> logaritmus je v tomto případě ideálním řešením.
hist(log(eurovolby$OBYVATEL), breaks = 50)

eurovolby$OBYVATEL_LOG <- log(eurovolby$OBYVATEL) # vytvoříme novou proměnnou

# Model nyní spočítáme s novou ("lepší") proměnnou.
# Regresní analýza:
reg_ucast_m3 <- lm(formula = UCAST ~ NEZAMESTNANYCH
                   + OBYVATEL_LOG + VEK + ZENY_PROC, data = eurovolby)
summary(reg_ucast_m3)

# Podívejme se na výsledky všech tří modelů:
stargazer(list(reg_ucast_m1, reg_ucast_m2, reg_ucast_m3), type = "html", out = "export/reg_ucast.html")

# Vysvětlovací schopnost modelu č. 3 se oproti předchozím zvýšila a
# hlavně jsme odhalili statistickou významnost počtu obyvatel.

# Provedeme simulaci:
range(eurovolby$OBYVATEL_LOG)

ci <- data.frame(NEZAMESTNANYCH = mean(eurovolby$NEZAMESTNANYCH), OBYVATEL_LOG = seq(2,15,0.1),
                 VEK = mean(eurovolby$VEK), ZENY_PROC = mean(eurovolby$ZENY_PROC))
# Následně odhalíme hodnoty závisle proměnné
sim <- predict(reg_ucast_m3, newdata = ci, interval = "confidence", level = 0.99)


plot(eurovolby$OBYVATEL_LOG, eurovolby$UCAST,
     ylab = "Volební účast (%)", xlab = "Počet obyvatel (log)",
     main = "Vztah velikosti obce a volební účasti (komplexní model)")
lines(seq(2,15,0.1),sim[,1], lt = "solid", lwd = 2, col = "orange")
lines(seq(2,15,0.1),sim[,2], lt = "dashed", lwd = 1, col = "orange")
lines(seq(2,15,0.1),sim[,3], lt = "dashed", lwd = 1, col = "orange")


# Pozor ale na interpretaci! Proměnná není ve své původní podobě, ale je zlogaritmovaná!

plot(eurovolby$OBYVATEL, eurovolby$UCAST,
     ylab = "Volební účast (%)", xlab = "Počet obyvatel",
     main = "Vztah velikosti obce a volební účasti (komplexní model)",
     xlim = c(0,10000))
lines(c(exp(seq(2,15,0.1))),sim[,1], lt = "solid", lwd = 2, col = "orange")
lines(c(exp(seq(2,15,0.1))),sim[,2], lt = "dashed", lwd = 1, col = "orange")
lines(c(exp(seq(2,15,0.1))),sim[,3], lt = "dashed", lwd = 1, col = "orange")

# Vidíme, že vztah díky transformaci už není lineární, ale právě logaritmický.



# Každou regresní analýzu by měly následovat analýzy robustnosti výsledků.
# Testuje se, jak moc jsou výsledky solidní - vytváříme různé modely s odlehlými hodnotami a bez nich,
# operacionalizujeme proměnné různým způsobem a sledujeme, zda vztah pořád přetrvává atd.




